library(haven)
library(psych)
library(plyr)
library(ggplot2)
library(cowplot)
library(Hmisc)
library(survey)
library(stargazer)
library(xtable)

setwd("~/data/lapop")
#setwd("c:/data/lapop")

# --- 2021 ---

# Load data
raw21 <- read_dta("lapop_rep/lapop_21_raw.dta")

# Strip variable and value labels
raw21 <- zap_label(raw21)
raw21 <- zap_labels(raw21)

# Create a pid variable that is 0 for other and independent and 1 for Dems, -1 for Reps
raw21$pid <- 0
raw21$pid[raw21$usvb1011 == 4001] <- -1  
raw21$pid[raw21$usvb1011 == 4002] <-  1  
table(raw21$usvb1011, raw21$pid)

# Data Collection June-July 2021
raw21$year <- 2021.6

# NOTE Variables used (unless noted before)
# jc13: Some people say that under some circumstances it would be justified for the military of this country to take
# power by a coup d’état (military coup). In your opinion, would a military coup be justified…
# When there is a lot of corruption.
# 1. Yes, it is justified / 2. No, it is not justified

# jc15a: Do you believe that when the country is facing very difficult times it is justifiable 
# for the president of the country to close the Congress and govern without Congress? 
# 1. Yes, it is justified / 2. No, it is not justified

#Add the variables missing
#raw21$jc13 <- NA
#raw21$jc15a <- NA
raw21$jc16a <- NA
raw21$e3 <- NA
raw21$pop101 <- NA
raw21$pop102 <- NA
raw21$pop103 <- NA

# define variables to keep
keeps21 <- c("jc13", "jc15a", "jc16a", "e3", "pop101", "pop102", "pop103", "pid", "wt", "year")

# Make new dataset with only kept variables
clean21 <- raw21[keeps21]

# Check new dataset
head(clean21)
psych::describe(clean21)

# --- 2019 ---

# See notes for 2021
raw19 <- read_dta("lapop_rep/lapop_19_raw.dta")
raw19 <- zap_label(raw19)
raw19 <- zap_labels(raw19)

raw19$pid <- 0
raw19$pid[raw19$usvb1011 == 4001] <- -1  
raw19$pid[raw19$usvb1011 == 4002] <-  1  
table(raw19$usvb1011, raw19$pid)

# Data Collection July 2019
raw19$year <- 2019.6

# NOTE new variables added as compared to prior datasets
# jc16a: Do you believe that when the country is facing very difficult times it is justifiable for the president/prime minister of the
# country to dissolve the Supreme Court of Justice and govern without the Supreme Court of Justice?
# 1. Yes, it is justified / 2. No, it is not justified

#Add the variables missing
raw19$jc13 <- NA
#raw19$jc15a <- NA
#raw19$jc16a <- NA
raw19$e3 <- NA
raw19$pop101 <- NA
raw19$pop102 <- NA
raw19$pop103 <- NA

# define variables to keep
keeps19 <- c("jc13", "jc15a", "jc16a", "e3", "pop101", "pop102", "pop103", "pid", "wt", "year")

clean19 <- raw19[keeps19]
psych::describe(clean19)

# --- 2017 ---

raw17 <- read_dta("lapop_rep/lapop_17_raw.dta")
raw17 <- zap_label(raw17)
raw17 <- zap_labels(raw17)

raw17$pid <- 0
raw17$pid[raw17$usvb1011 == 4001] <- -1  
raw17$pid[raw17$usvb1011 == 4002] <-  1  
table(raw17$usvb1011, raw17$pid)

# Data Collection mid-May 2017
raw17$year <- 2017.5

#Add the variables missing
#raw17$jc13 <- NA
#raw17$jc15a <- NA
raw17$jc16a <- NA
raw17$e3 <- NA
raw17$pop101 <- NA
raw17$pop102 <- NA
raw17$pop103 <- NA

# define variables to keep
keeps17 <- c("jc13", "jc15a", "jc16a", "e3", "pop101", "pop102", "pop103", "pid", "wt", "year")

clean17 <- raw17[keeps17]

head(clean17)
psych::describe(clean17)

# --- 2014 ---

raw14 <- read_dta("lapop_rep/lapop_14_raw.dta")
raw14 <- zap_label(raw14)
raw14 <- zap_labels(raw14)

raw14$pid <- 0
raw14$pid[raw14$usvb1011 == 4001] <- -1  
raw14$pid[raw14$usvb1011 == 4002] <-  1  
table(raw14$usvb1011, raw14$pid)

# Data Collection end of June - early July 2014
raw14$year <- 2014.6

# NOTE new variables added as compared to prior datasets
# e3:  Of people participating in a group working to violently overthrow an elected government.
# How much do you approve or disapprove? 
# 1 strongly disapprove to 10 strongly approve 

#Add the variables missing
#raw14$jc13 <- NA
#raw14$jc15a <- NA
raw14$jc16a <- NA
#raw14$e3 <- NA
raw14$pop101 <- NA
raw14$pop102 <- NA
raw14$pop103 <- NA

# define variables to keep
keeps14 <- c("jc13", "jc15a", "jc16a", "e3", "pop101", "pop102", "pop103", "pid", "wt", "year")

clean14 <- raw14[keeps14]

head(clean14)
psych::describe(clean14)

# --- 2012 ---

raw12 <- read_dta("lapop_rep/lapop_12_raw.dta")
raw12 <- zap_label(raw12)
raw12 <- zap_labels(raw12)

# Note variable name change for original party id
raw12$pid <- 0
raw12$pid[raw12$vb11s == 2] <- -1  
raw12$pid[raw12$vb11s == 1] <-  1  
table(raw12$vb11s, raw12$pid)

# Data Collection April 2012
raw12$year <- 2012.3

# NOTE new variables added as compared to prior datasets
# pop101: . It is necessary for the progress of this country that our presidents limit the voice
# and vote of opposition parties, how much do you agree or disagree with that view?
# 1 strongly disapprove to 10 strongly approve 

#Add the variables missing
#raw12$jc13 <- NA
#raw12$jc15a <- NA
#raw12$jc16a <- NA
#raw12$e3 <- NA
#raw12$pop101 <- NA
raw12$pop102 <- NA
raw12$pop103 <- NA

# define variables to keep
keeps12 <- c("jc13", "jc15a", "jc16a", "e3", "pop101", "pop102", "pop103", "pid", "wt", "year")

clean12 <- raw12[keeps12]
psych::describe(clean12)

# --- 2010 ---

raw10 <- read_dta("lapop_rep/lapop_10_raw.dta")
raw10 <- zap_label(raw10)
raw10 <- zap_labels(raw10)

raw10$pid <- 0
raw10$pid[raw10$vb11s == 2] <- -1  
raw10$pid[raw10$vb11s == 1] <-  1  
table(raw10$vb11s, raw10$pid)

# Data Collection second half of March in 2010 
raw10$year <- 2010.3

summary(raw10$wt)

# NOTE new variables added as compared to prior datasets
# pop102: When the Congress hinders the work of our government, our presidents should
# govern without the Congress. How much do you agree or disagree with that view? 
# 1 strongly disapprove to 10 strongly approve 

# pop103: When the Supreme Court blocks the work of our government, the Supreme Court
# should be disregarded by our presidents. How much do you agree or disagree with that view?
# 1 strongly disapprove to 10 strongly approve 

# define variables to keep (this  year is the one that has it all)
keeps10 <- c("jc13", "jc15a", "jc16a", "e3", "pop101", "pop102", "pop103", "pid", "wt", "year")

clean10 <- raw10[keeps10]
psych::describe(clean10)

# --- 2008 ---

raw08 <- read_dta("lapop_rep/lapop_08_raw.dta")
raw08 <- zap_label(raw08)
raw08 <- zap_labels(raw08)

# Note variable name change for original party id
raw08$pid <- 0
raw08$pid[raw08$pid3 == 2] <- -1  
raw08$pid[raw08$pid3 == 1] <-  1  
table(raw08$pid3, raw08$pid)

# Data Collection April 2008
raw08$year <- 2008.3

#Add the variables missing
raw08$jc13 <- NA
raw08$jc15a <- NA
raw08$jc16a <- NA
#raw08$e3 <- NA
#raw08$pop101 <- NA
#raw08$pop102 <- NA
#raw08$pop103 <- NA

# define variables to keep
keeps08 <- c("jc13", "jc15a", "jc16a", "e3", "pop101", "pop102", "pop103", "pid", "wt", "year")

clean08 <- raw08[keeps08]
psych::describe(clean08)

# --- 2006 ---

raw06 <- read_dta("lapop_rep/lapop_06_raw.dta")
raw06 <- zap_label(raw06)
raw06 <- zap_labels(raw06)

# Note variable name change for original party id (caps)
raw06$pid <- 0
raw06$pid[raw06$VB11 == 2] <- -1  
raw06$pid[raw06$VB11 == 1] <-  1  
table(raw06$VB11, raw06$pid)

# No weights supplied. Each individual has the same weight.
raw06$wt <- 1

# Data Collection took place God only knows when. Approximating with mid-year. 
raw06$year <- 2006.5
  
# need to redo variable names to match other years (caps and add "a" at the end)
raw06$jc15a <- raw06$JC15
raw06$jc16a <- raw06$JC16

#Add the variables missing
raw06$jc13 <- NA
#raw06$jc15a <- NA
#raw06$jc16a <- NA
raw06$e3 <- NA
raw06$pop101 <- NA
raw06$pop102 <- NA
raw06$pop103 <- NA

# define variables to keep
keeps06 <- c("jc13", "jc15a", "jc16a", "e3", "pop101", "pop102", "pop103", "pid", "wt", "year")

clean06 <- raw06[keeps06]
psych::describe(clean06)

dhlapop <- rbind(clean06, clean08, clean10, clean12, clean14, clean17, clean19, clean21)

head(dhlapop)
psych::describe(dhlapop)

# --- Recode Work ---

# Recode jc13 to go from 0 to 100
dhlapop$coupcorrupt <- (2 - dhlapop$jc13) * 100
table(dhlapop$coupcorrupt, dhlapop$jc13)

# Recode jc15a to go from 0 to 100
dhlapop$dissolve <- (2 - dhlapop$jc15a) * 100
table(dhlapop$dissolve, dhlapop$jc15a)

# Recode jc16a to go from 0 to 100
dhlapop$supcor <- (2 - dhlapop$jc16a) * 100
table(dhlapop$supcor, dhlapop$jc16a)

# Recode e3 to go from 0 (less democracy eroding) to 100 (more democracy eroding)
dhlapop$violentoverthrow <- (dhlapop$e3 - 1) / 9 * 100
table(dhlapop$violentoverthrow, dhlapop$e3)

# Recode pop101, 102, 103 to go from 0 to 100
dhlapop$limitopposition <- (dhlapop$pop101 - 1) / 6 * 100
table(dhlapop$limitopposition, dhlapop$pop101)

dhlapop$limitcongress <- (dhlapop$pop102 - 1) / 6 * 100
table(dhlapop$limitcongress, dhlapop$pop102)

dhlapop$limitsupcor <- (dhlapop$pop103 - 1) / 6 * 100
table(dhlapop$limitsupcor, dhlapop$pop103)

psych::describe(dhlapop)

## Unweighted
## Create data frame for plotting
#lapops <- ddply(dhlapop, .(pid, year), summarise,
#                discongress.mean      = mean(dissolve, na.rm = T),
#                discongress.se        = sd(dissolve, na.rm = T) / sqrt(length(dissolve)),
#                dissupcor.mean        = mean(supcor, na.rm = T),
#                dissupcor.se          = sd(supcor, na.rm = T) / sqrt(length(supcor)),
#                coupcorrupt.mean      = mean(coupcorrupt, na.rm = T),
#                coupcorrupt.se        = sd(coupcorrupt, na.rm = T) / sqrt(length(coupcorrupt)),
#                violentoverthrow.mean = mean(violentoverthrow, na.rm = T),
#                violentoverthrow.se   = sd(violentoverthrow, na.rm = T) / sqrt(length(violentoverthrow)),
#                limitopposition.mean  = mean(limitopposition, na.rm = T),
#                limitopposition.se    = sd(limitopposition, na.rm = T) / sqrt(length(limitopposition)),
#                limitcongress.mean    = mean(limitcongress, na.rm = T),
#                limitcongress.se      = sd(limitcongress, na.rm = T) / sqrt(length(limitcongress)),
#                limitsupcor.mean      = mean(limitsupcor, na.rm = T),
#                limitsupcor.se        = sd(limitsupcor, na.rm = T) / sqrt(length(limitsupcor))
#                )
#lapops <- subset(lapops, pid %in% c(-1, 1))

# Weighted
# Create data frame for plotting
lapops <- ddply(dhlapop, .(pid, year), summarise,
                discongress.mean      = weighted.mean(dissolve, wt, na.rm = T),
                discongress.se        = sqrt(wtd.var(dissolve, wt, na.rm = T) / length(na.omit(dissolve))),
                dissupcor.mean        = weighted.mean(supcor, wt, na.rm = T),
                dissupcor.se          = sqrt(wtd.var(supcor, wt, na.rm = T) / length(na.omit(supcor))),
                coupcorrupt.mean      = weighted.mean(coupcorrupt, wt, na.rm = T),
                coupcorrupt.se        = sqrt(wtd.var(coupcorrupt, wt, na.rm = T) / length(na.omit(coupcorrupt))),
                violentoverthrow.mean = weighted.mean(violentoverthrow, wt, na.rm = T),
                violentoverthrow.se   = sqrt(wtd.var(violentoverthrow, wt, na.rm = T) / length(na.omit(violentoverthrow))),
                limitopposition.mean  = weighted.mean(limitopposition, wt, na.rm = T),
                limitopposition.se    = sqrt(wtd.var(limitopposition, wt, na.rm = T) / length(na.omit(limitopposition))),
                limitcongress.mean    = weighted.mean(limitcongress, wt, na.rm = T),
                limitcongress.se      = sqrt(wtd.var(limitcongress, wt, na.rm = T) / length(na.omit(limitcongress))),
                limitsupcor.mean      = weighted.mean(limitsupcor, wt, na.rm = T),
                limitsupcor.se        = sqrt(wtd.var(limitsupcor, wt, na.rm = T) / (length(na.omit(limitsupcor))))
)
lapops <- subset(lapops, pid %in% c(-1, 1))

weighted.mean(dhlapop$dissolve, dhlapop$wt, na.rm = T)
mean(dhlapop$dissolve, na.rm = T)

# Make factors from party variable
lapops$party<-factor(lapops$pid,
                     levels = c(-1,1),
                     labels = c( "R","D"))

# This is for the shading of the circle around the letter in the plot
# probably could make it the background sharing as well. Would be more elegant. But I can't be bothered now.
lapops$fill <- "a"
lapops$fill[lapops$year > 2007 & lapops$year < 2009] <- "b"
lapops$fill[lapops$year > 2009 & lapops$year < 2011] <- "c"
lapops$fill[lapops$year > 2011 & lapops$year < 2017] <- "d"
lapops$fill[lapops$year > 2017 & lapops$year < 2019] <- "a"
lapops$fill[lapops$year > 2019 & lapops$year < 2021] <- "b"
lapops$fill[lapops$year > 2021] <- "c"

write.csv(lapops, file = "plotdata.csv")

lapops

# OK. To make life easy, I will make a plotdata for each plot. These will remove the rows with NAs for the variable we are 
# plotting. Otherwise the line disappears when there is an NA in between two data points (and it's a nightmare to compensate)
plotdata <- lapops[complete.cases(lapops$discongress.mean), ]
pd <- position_dodge(0.2)

series.legislature <- ggplot(plotdata, aes(x = year, y = discongress.mean, group = party, color = party, label = party)) + 
  geom_rect(aes(xmin = -Inf, xmax = 2007, ymin = -Inf, ymax = 0),   fill = "grey80", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = -Inf, xmax = 2007, ymin = 0,    ymax = Inf), fill = "grey92", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2007, xmax = 2009, ymin = 0,    ymax = Inf), fill = "grey92", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2007, xmax = 2009, ymin = -Inf, ymax = 0),   fill = "grey65",   color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2009, xmax = 2011, ymin = -Inf, ymax = 0),   fill = "grey80", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2009, xmax = 2011, ymin = 0,    ymax = Inf), fill = "white", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2011, xmax = 2017, ymin = -Inf, ymax = 0),   fill = "grey65",   color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2011, xmax = 2017, ymin = 0,    ymax = Inf), fill = "white", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2017, xmax = 2019, ymin = -Inf, ymax = 0),   fill = "grey80", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2017, xmax = 2019, ymin = 0,    ymax = Inf), fill = "grey92", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2019, xmax = 2021, ymin = -Inf, ymax = 0),   fill = "grey65", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2019, xmax = 2021, ymin = 0,    ymax = Inf), fill = "grey92", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2021, xmax = Inf,  ymin = -Inf, ymax = 0),   fill = "grey80", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2021, xmax = Inf,  ymin = 0,    ymax = Inf), fill = "white", color = NA, alpha = 0.3) +
  geom_errorbar(aes(ymin = discongress.mean - 2 * discongress.se, 
                    ymax = discongress.mean + 2 * discongress.se), 
                width = 0, position = pd) +
  geom_line(linewidth = 0.2) +
  geom_label(position = pd, size = 2.5, label.r = unit(.45, "lines"), label.size = 0, aes(fill = fill)) +
  scale_fill_manual(values = c("grey92", "grey92", "white", "white")) +
  scale_x_continuous(limits = c(2006, 2022), breaks = seq(2006, 2022, by = 2)) +
  scale_y_continuous(limits = c(0, 100)) + # was 57
  ylab("Sometimes justifiable for the president to close the Congress") + xlab("Year") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = "none", strip.background = element_rect(fill = "white")) +
  scale_color_manual(name = "Party",
                     values = c('grey10', 'grey45'),
                     labels = c("Republicans", "Democrats")) 

print(series.legislature)

ggsave("timeseriesFINALRR.png", plot = series.legislature, width = 5, height = 5.7, limitsize = F, dpi = 1200)
ggsave("timeseriesFINALRR.tiff", plot = series.legislature, width = 5, height = 5.7, limitsize = F, dpi = 300)
ggsave("timeseriesFINALRR.pdf", plot = series.legislature, width = 5, height = 5.1)
ggsave("timeseriesFINALRR.eps", plot = series.legislature, width = 5, height = 5.7, limitsize = F, dpi = 1200)

plotdata <- lapops[complete.cases(lapops$dissupcor.mean), ]
print(plotdata)

series.supcor <- ggplot(plotdata, aes(x = year, y = dissupcor.mean, group = party, color = party, label = party)) + 
  geom_rect(aes(xmin = -Inf, xmax = 2009, ymin = -Inf, ymax = Inf), fill = "grey92",  color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2009, xmax = 2017, ymin = -Inf, ymax = Inf), fill = "white", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2017, xmax = 2021, ymin = -Inf, ymax = Inf), fill = "grey92",  color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2021, xmax = Inf,  ymin = -Inf, ymax = Inf), fill = "white", color = NA, alpha = 0.3) +
  geom_errorbar(aes(ymin = dissupcor.mean - 2 * dissupcor.se, 
                    ymax = dissupcor.mean + 2 * dissupcor.se), 
                width = 0, position = pd) +
  geom_line(linewidth = 0.2) +
  geom_label(position = pd, size = 2.5, label.r = unit(.45, "lines"), label.size = 0, aes(fill = fill)) +
  scale_fill_manual(values = c("grey92", "grey92", "white", "white")) +
  scale_x_continuous(limits = c(2006, 2022), breaks = seq(2006, 2022, by = 2)) +
  scale_y_continuous(limits = c(0, 100)) +
  ylab("Sometimes justifiable for the president to close the Supreme Court") + xlab("Year") +  
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = "none", strip.background = element_rect(fill = "white")) +
  scale_color_manual(name = "Party",
                     values = c('grey10', 'grey45'),
                     labels = c("Republicans", "Democrats"))

print(series.supcor)

ggsave("AppendixC1RR.png", plot = series.supcor, width = 5, height = 5.7, limitsize = F, dpi = 1200)
ggsave("AppendixC1RR.tiff", plot = series.supcor, width = 5, height = 5.7, limitsize = F, dpi = 300)
ggsave("AppendixC1RR.pdf", plot = series.supcor, width = 5, height = 5.1)

plotdata <- lapops[complete.cases(lapops$limitopposition.mean), ]

series.limitopposition2 <- ggplot(plotdata, aes(x = year, y = limitopposition.mean, group = party, color = party, label = party)) + 
  geom_rect(aes(xmin = -Inf, xmax = 2009, ymin = -Inf, ymax = Inf), fill = "grey92",  color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2009, xmax = 2017, ymin = -Inf, ymax = Inf), fill = "white", color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2017, xmax = 2021, ymin = -Inf, ymax = Inf), fill = "grey92",  color = NA, alpha = 0.3) +
  geom_rect(aes(xmin = 2021, xmax = Inf,  ymin = -Inf, ymax = Inf), fill = "white", color = NA, alpha = 0.3) +
  geom_errorbar(aes(ymin = limitopposition.mean - 2 * limitopposition.se, 
                    ymax = limitopposition.mean + 2 * limitopposition.se), 
                width = 0, position = pd) +
  geom_line(size = 0.2) +
  geom_label(position = pd, size = 2.5, label.r = unit(.45, "lines"), label.size = 0, aes(fill = fill)) +
  scale_fill_manual(values = c("grey92", "white", "white")) +
  scale_x_continuous(limits = c(2006, 2022), breaks = seq(2006, 2022, by = 2)) +
  scale_y_continuous(limits = c(0, 100)) +
  ylab("For progress Presidents can limit the voice & vote of the opposition") + xlab("Year") +  
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        legend.position = "none", strip.background = element_rect(fill = "white")) +
  scale_color_manual(name = "Party",
                     values = c('grey10', 'grey45'),
                     labels = c("Republicans", "Democrats"))

print(series.limitopposition2)

ggsave("AppendixC2RR.png", plot = series.limitopposition2, width = 5, height = 5.7, limitsize = F, dpi = 1200)
ggsave("AppendixC2RR.tiff", plot = series.limitopposition2, width = 5, height = 5.7, limitsize = F, dpi = 300)
ggsave("AppendixC2RR.pdf", plot = series.limitopposition2, width = 5, height = 5.1)

lapops
(tableInAp <- cbind(lapops[, 2:6], lapops[, 11:12], lapops[, 17]))
tableInAp

print(xtable(tableInAp), include.rownames = FALSE)

psych::describe(dhlapop)

table(dhlapop$year)

dhlapop$respr <- 0
dhlapop$respr[dhlapop$pid == -1] <- 1
dhlapop$respi <- 0
dhlapop$respi[dhlapop$pid == 0]  <- 1
dhlapop$respd <- 0
dhlapop$respd[dhlapop$pid == 1]  <- 1
dhlapop$presd <- 1
dhlapop$presd[dhlapop$year == 2006.5] <- 0
dhlapop$presd[dhlapop$year > 2016] <- 0
dhlapop$presd[dhlapop$year > 2020] <- 1
dhlapop$presr <- 1 - dhlapop$presd
dhlapop$divg <- 0
dhlapop$divg[dhlapop$year == 2008.3] <- 1
dhlapop$divg[dhlapop$year == 2012.3] <- 1
dhlapop$divg[dhlapop$year == 2014.6] <- 1
dhlapop$divg[dhlapop$year == 2019.6] <- 1

dhlapop$yearid <- as.numeric(as.factor(dhlapop$year))
dhlapop$yearid <- as.factor(dhlapop$yearid)

table(dhlapop$year, dhlapop$presd)
table(dhlapop$year, dhlapop$divg)

dhlapop2 <- dhlapop[dhlapop$respi == 0, ]
table(dhlapop2$yearid, dhlapop2$divg)

dhlapop2$same <- ifelse(dhlapop2$presd == dhlapop2$respd, 1, 0)

regdata3 <- svydesign(ids = ~yearid, weights = ~wt, data = dhlapop2) 
dh3.res1 <- svyglm(dissolve ~ same, design = regdata3, family = gaussian(link = "identity"))
dh3.res2 <- svyglm(dissolve ~ same * divg, design = regdata3, family = gaussian(link = "identity"))
stargazer(dh3.res1, dh3.res2, type = "text", report = ('vc*p'))

stargazer(dh3.res1, dh3.res2, report = ('vc*p'))
